Normalization integrals of orthogonal Heun functions 



Peter A. Becker 



Center for Earth Observing and Space Research 
Institute for Computational Sciences and Informatics, 
and Department of Physics and Astronomy, 
George Mason University, Fairfax, VA 22030-4444 



(pbeckerOgmu . edu) 

A formula for evaluating the quadratic normalization integrals of orthogonal Heun 
functions over the real interval < a; < 1 is derived using a simple limiting procedure 
based upon the associated differential equation. The resulting expression gives the value 
of the normalization integral explicitly in terms of the local power-series solutions about 
X = and x = 1 and their derivatives. This provides an extremely efficient alternative to 
nimierical integration for the development of an orthonormal basis using Heun functions, 
because all of the required information is available as a by-product of the search for the 
eigenvalues of the differential equation. 

PACS numbers: 02.30.Gp, 02.30.Hq, 02.70.-c, 02.60.Jh 

Byline: Normalization integrals of Heun functions 

(Accepted for publication in the Journal of Mathematical Physics) 



1 



I. INTRODUCTION 

Heun's equation^ is the most general Fuchsian equation of second order with four 
regular singular points, and it is therefore of considerable importance in mathematical 
physics. Special cases of the Heun equation include the hypergeometric, confluent hyper- 
geometric, Lame, Bessel, Legendre, and Laguerre equations. As a practical matter, the 
most important solutions to the Heun equation are those orthogonal functions satisfying 
prescribed boundary conditions at two adjacent singular points.^ The development of an 
orthonormal basis using these functions is a two-step process. The first step is to search 
for the eigenvalues, and the second is to evaluate the quadratic normalization integrals of 
the associated eigenfunctions, which are orthogonal Heun functions. The normalization 
integrals are usually evaluated numerically, which is not an especially efficient procedure 
given the nature of the eigenfunctions and the associated series representations. In order to 
provide a useful alternative to numerical integration, in this paper we derive a new formula 
for directly evaluating the quadratic normalization integrals of orthogonal Heun functions 
over the real interval [0, 1]. The formula obtained (eq. 29) utilizes only information avail- 
able as a by-product of the search for the eigenvalues, and therefore greatly improves both 
the efficiency and the accuracy of numerical procedures involving Heun functions. 

II. HEUN FUNCTIONS 

We begin by writing Heun's equation in the standard form first adopted by Erdelyi et 

al.,3 

dx"^ \x X — 1 X — a J dx x{x — l){x — a)^ ' 
represented by the Riemann P-symbol 

f 1 a oo ] 
pI a x \ , (2) 
[l-7 1-S 1-e P J 

with regular singular points located at x = 0, 1, a, oo. The parameter A falls outside the 

domain of the usual Riemann classification scheme, and is therefore referred to as an acces- 
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sory (or aiixiliary) parameter. In many applications, A plays the role of an eigenparameter. 
The five exponent parameters a, /3, 7, S, e are connected via Riemann's relation 

a + /?-7-5-e + l = 0, (3) 

and therefore only four of them are independent. The total number of free parameters is 
six, and this number cannot be reduced by any transformation. 

Heun^ used the method of Frobenius to derive local powcr-scrics solutions to (1), 
generating a three-term recursion relation for the expansion coefficients. Two linearly 
independent power-series (Frobenius) solutions exist in the neighborhood of any one of the 
singular points, and, in general, analytic continuation of a single Frobenius solution about 
one singularity into the neighborhood of a second, adjacent singularity generates a linear 
combination of the two local Frobenius solutions about the second singularity. 

The most important solutions are those that are simultaneously local Frobenius so- 
lutions about two adjacent singular points. These are referred to as Heun functions, and 
often arise when physical boundary conditions are applied to solutions of the differential 
equation. Taking the parameters a, /3, 7, 5, a to be constants, the problem of finding 
a Heun function becomes a singular Sturm-Liouville eigenvalue problem for A. However, 
since no formula is available for performing the analytic continuation between two adjacent 
singularities in the case of the general Heun equation (in contrast to the subcase of the 
hypergeometric equation), no general closed-form expression for the eigenvalues A^ exists. 
Heun functions are infinite series in general, although in special cases the series truncates, 
leaving a Heun polynomial. 

There are four classes of Heun functions for a given pair of adjacent singularities 
X = So and x = si that are distinguished by the values of the corresponding exponents 
((To, CTi). In this paper, we set so = and si = 1, and classify the Heun functions according 
to the usual scheme based upon the values of the associated exponents: class I (0, 0); class 
II (1 — 7, 0); class III (0, 1 — 5); and class IV (1 — 7, 1 — 5). We focus here on the behavior 

3 



of Heun functions within the real interval < a; < 1, and we assume throughout that 
a^[0, 1]. 

III. DETERMINATION OF THE EIGENVALUES 

Let the function yo{X,x) be a local Frobenius solution of (1) in the neighborhood of 
X = 0, and let the function j/i(A, x) be a local Frobenius solution in the neighborhood of 
X = 1. The properties of these solutions, including the recursion relation for the expansion 
coefficients, have been fully discussed.^ Let us suppose that these functions are normalized 
according to the prescription adopted by Heun, so that 

lim x-'^" yo{X, x) = lim {x - yi{X, x) ^ 1 . (4) 

a;— »0 x—^l 

It is known from the basic theory of these solutions that in general yo converges for < 
min(l, \a\) and yi converges for |a; — 1| < min(l, |a — 1|). Hence when a ^ [0, 1] as assumed 
here, there exists a region of mutual convergence, within which both yo ^ind j/i converge. 
The region of mutual convergence is contained within the interval [0, 1]. 
The Wronskian of yo and yi is defined by 

W{X,x) = yo{X,x)^{X,x) -yi{X,x)^{X,x) . (5) 

In order to obtain a Heun function, the Wronskian must vanish, and therefore the eigen- 
value equation for A„ becomes 

W{Xn,x)=0. (6) 

When this condition is satisfied, yQ{XnjX) and yi{Xn,x) are linearly dependent functions 
(although not equal in general), and the solution to (1) with A = is the Heun func- 
tion Hn{x). We set the normalization of Hn{x) by stipulating that Hn{x) = yo{Xn,x) 
in the neighborhood of a: = 0. In the neighborhood of a; = 1, we find that Hn{x) = 
A{Xn) yi{Xn, x), where the value of A{Xn) is determined by requiring that Hn{x) be con- 
tinuous at an arbitrary point within the region of mutual convergence of yo and yi. 
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We can establish the functional form of the Wronskian by examining the self-adjoint 
version of the Heun equation, 

Cy-Xuj{x)y = Q, (7) 
where the weight function u3{x) is defined by 

uj{x) = x^-^{x-lY-^{x-ay-^, (8) 

and the operator C is defined by 

fill 1 

-\- a (5 xijo(x)y . (9) 
Since yo and yi are each solutions of (1) for the same value of A, we may write 

yo [C - \uj{x)] yi - yi [C - \uj{x)] yo = . (10) 

This yields an equation for the Wronskian; 

1 dW 'y 5 e 
+ T + 



Cy^- 
dx 



x^ix — V)(x — (x) 

dx 



W dx X X — 1 X — a' 
with solution 

W{X, x) - D{X) x-'^ix - iy\x -a)-', (11) 

where -D(A) is an unknown function. Note that the Wronskian vanishes for D{Xn) = 0, 
which is independent of x. Hence when we evaluate W{X, x) using (5) in order to calculate 
the eigenvalues using (6), we are free to pick any convenient value for x that lies within 
the region of mutual convergence of j/o ^ind j/i. 

IV. ORTHOGONALITY RELATIONS 

Let Hn{x) and Hm{x) be Heun functions of the same class associated with eigenvalues 
Xn and A^. Since they are each solutions of (1), we have, by analogy with (10), 

Hn [C - Xm i^{x)] Hm -Hm[C- A„ Uj{x)] Hn = . (12) 
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After simplifying and integrating over the interval [0, 1], we obtain 



dx dx 



1 







(13) 



where 

p{x) = x^{x-lf{x-aY. (14) 

Based upon the asymptotic behavior of the local Probenius solutions j/o ^ind j/i, we find 
that the right-hand side of (13) vanishes when one of the following sets of class-dependent 
conditions is satisfied; 

class I: 3?7 > , 3f?5 > 

class II: 3?7<2, 3?(5>0 , . 

class III: 3fJ7>0, 3?5<2 ^ ' 

class IV: 3^7 < 2 , m<2. 

We shall refer to (15) as the set of existence conditions for orthogonal Heun functions on 

the interval [0, 1]. When these conditions are met, we obtain the standard orthogonality 

relation 

(An - Am) / ijo{x) Hn{x) Hm{x) dx = Q . (16) 

V. NORMALIZATION INTEGRALS 

Naturally, the integral in (16) does not vanish when n = m, and in this case it is 
necessary to establish its value in order to develop an orthonormal basis using orthogonal 
Heun functions. Numerical integration is always available as an option, but this is a 
very inefficient approach to the problem. Other procedures have been devised, the most 
interesting being the method developed by Erdelyi,^ which is based upon expansions of 
Heun functions as series of degenerate hypergeometric functions (Jacobi polynomials). 
Lambe and Ward^ developed a technique for evaluating normalization integrals for Heun 
polynomials, but these results are not applicable to the more general (and much more 
common) case of Heun functions. In this section we develop a new formula for the explicit 
evaluation of these integrals. 
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We can derive a formula for evaluating the quadratic normalization integral 



In= ^(x) [Hn{x)f dx 

Jo 



(17) 



by generalizing the approach taken in the preceding section. Proceeding as before, we note 
that since yo{X,x) and yo{Xn,x) are each Frobenius solutions of (1) in the neighborhood 
of x = 0, it follows by analogy with (12) that 



yoiX, x) [C - Xn uj{x)] yo{Xn, x) - yo{Xn, x) [C - Xuj{x)] yo{X, x) = , 



(18) 



where An is an eigenvalue and A is arbitrary. After simplifying and integrating with respect 
to X, we now obtain 

(A -An) / oj{x')yo{X,x')yo(Xn,x')dx' 
Jo 

= p{x) 



yo{Xn,x) -^\\^) - yo{X,x) -g^y^n,x) 



(19) 



where the right-hand side vanishes as a; — > provided the appropriate set of existence 
conditions in (15) is satisfied. 

Upon examination of (19), we find that both sides of the equation vanish as A — > An- 
We can therefore establish the value of the indefinite integral in the limit A ^ A^ using 
L'Hopital's rule, which yields 

uj{x') [yo{Xn,x')]'^ dx' 



I 

Jo 



= p{x) 



(20) 



To obtain a formula for the desired normalization integral In, we must let a; ^ 1 in (20), 
and this requires analytic continuation of yo into the neighborhood of a; = 1. For general 
values of A, we write the analytic continuation of yo as 



yo(A, x) = A{X) yi{X, x) + B{X) yi{X, x) , 



(21) 



where yi is the Frobenius solution about a; = 1 with exponent zero, and yi is the Frobenius 
solution about a; = 1 with exponent 1 — S. At this point we shall restrict our attention to 



Heun functions of classes I or II, so that the exponent at a; = 1 is zero. This restriction 
will be removed later. Note that B{Xn) must vanish so that we obtain a class I or II Heun 
function when A = A^. Substituting into (20), we obtain after some algebra 



In= Uj{x) [Hn{x)f dx 

Jo 



= lim 



2.. .r. X / .r. xn (22) 



. / . d^yi dB dyi \ . dyi ( . dyi dB ^ 

^y^^'^dxd^+dx-d^-'^^r^+dxy' 



p 



where we have used the fact that in the neighborhood of x = 1, Hn{x) is the analytic 
continuation of yo{Xn, x), and in the neighborhood of a; = 0, Hn{x) is identical to yo{Xni x). 

Based upon asymptotic analysis of yi and yi, we find that in the limit a; —> 1, all of 
the terms on the right-hand side of (22) vanish except the second one, so that we are left 
with 

In = lim p(a;)^(An) -7T-(An) ^— (An,a;)?/i(An,a;) . (23) 

x^i a A ox 

Hence we need only evaluate ^(An) and ^(A^) in terms of known functions in order to 
obtain our final result for the normalization integral /„. Evaluation of ^(An) is a simple 
matter, since the continuity of Hn{x) requires that 

A(A.) = (24) 

for any x within the region of mutual convergence of yo and yi. Evaluation of ^(A„) is 
slightly more complicated. We begin by differentiating (21) with respect to x to obtain 

^(A, x) = A{X) ^(A, x) + B{X) ^(A, x) . (25) 

Solving (21) and (25) for B{X) yields 



BiX) = . (26) 

^1 dx dx 



Bearing in mind that -B(A„) = 0, we find upon differentiating (26) with respect to A that 



where the Wronskian W{X,x) is defined by (5). We can use (24) and (27) respectively to 

dB 
dX 



replace A(A„) and 4x(An) in (23), yielding 



In — lini p— ^ g — yi 



dW 

I /7I-1 

(28) 

A — An 



The most singular term in the denominator is the one containing , and therefore in the 
limit a; — > 1 our final result for the normalization integral becomes 

l\{x) \lin{.x)f dx = -p{x) ^(An,x) y^i^ . (29) 

In passing to this expression, we have made use of the fact that ^^{Xn,x)p{x) and 
yo{Xnix)/yi{Xn,x) are both independent of x. Hence the right-hand side of (29) is ac- 
tually an invariant, which may be evaluated at any point within the region of mutual 
convergence of the local power-series solutions yo and yi. 

VI. EXTENSION TO HEUN FUNCTIONS OF CLASSES III AND IV 

We have derived our main result (29) under the assumption that Hn{x) is a class I or 
II Heun function. In the case of a Heun function of class III or IV, the exponent at a; = 1 
is 1 — 5, and the derivation presented above must be modified slightly. First we note that 
(20) remains valid because it is written in terms of the local solution y^ about a: = 0. The 
analytic continuation of j/o into the neighborhood of a; = 1 is still performed using (21), 
except now we interchange the definitions of j/i and yi, so that yi is the local Frobenius 
solution about a; = 1 with exponent 1 — 5 and yi is the local solution about a; = 1 with 
exponent zero. In this case we again require -B(A^) = 0, and following the same procedure 
as before we regain expression (22) . Due to the fact that the definitions of j/i and yi have 
been interchanged, in the limit a; — > 1, we find by asymptotic analysis that only the final 
term in (22) now contributes, so that in this case we obtain 

dB dyi 

= lim -p{x) A{Xn) —{Xn) ^{Xn,x)yi{Xn,x) . (30) 

a;-»l dX ox 
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Evaluation of ^(A„) and 4f (An) proceeds exactly as before, and we simply regain expres- 



sions (24) and (27). Using these results in (30) now yields 

' • (31) 



In = -P — , g^^ g~- yi 



A — A„ 

This is similar to (28), except now the most singular term in the denominator is the one 
containing and therefore in the limit x — > 1 we obtain 

In = -p{x) ^^{K,x) ^"^^" "' ^ , 

which is identical to (29) . We are therefore led to the following general conclusion: Equation 
(29) holds for Heun functions of any class, provided the existence conditions (15) are 
satisfied. We discuss the significance of this result and its natural role in computational 
algorithms below. 

VII. COMPUTATIONAL IMPLICATIONS 

Equation (29) has important practical consequences for the design of numerical al- 
gorithms used to develop orthonormal systems based on Heun functions, which are the 
solutions of interest in many mathematical and physical situations. In such cases the first 
step towards solution is the search for the associated eigenvalues A^- This search must 
proceed numerically in general, due to the lack of a closed-form expression for A„ in terms 
of the parameters a, 7, 5, a. The eigenvalues are obtained by isolating the roots of 
the Wronskian in (6). The root finding usually proceeds via Newton's method or possibly 
some more sophisticated algorithm. Most of these techniques require the evaluation of W 
and in order to generate a revised estimate of the true root A„, and the evaluation of 
these functions in turn involves the determination of the quantities 

dyo dyo d'^yo dyi dyi d'^yi 

yo ) Q ' a\ 1 ci\ ci 1 yi 1 Q. 1 Q,\ 1 Qi\ c 1 
ox oX oXox OX OA OXOX 

at each iteration. The values of j/o and J/i can be obtained using the well-known power- 
series representations, and the values of the derivatives can be obtained using term- by-term 
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differentiation. It is straiglitforward to demonstrate tliat tlie radii of convergence of the 
series for the derivatives are identical to those for the corresponding fundamental series, 
which are discussed in § II. 

Once the eigenvalues have been determined to acceptable precision, one generally 
needs to evaluate the associated quadratic normalization integrals in order to develop 
a set of orthonormal basis functions hn{x) using 



The conventional approach to the problem of evaluating is to integrate (17) numerically, 
in which case thousands of evaluations of Hn{x) would generally be required in order to 
establish the value of In to reasonable accuracy. However, such an inefficient procedure is 
no longer necessary with the availability of (29), because it allows the determination of to 
high precision using only the values of yo, j/i, and obtained in the final iteration of the 
root- finding stage of the algorithm. Hence no substantial additional calculation is necessary 
in order to determine In- The computational time required for developing orthonormal 
systems of Heun functions can therefore be reduced by several orders of magnitude by 
using (29) instead of numerical integration. 

We close by making a comparison between the method for evaluating In outlined here 
and that suggested by Erdelyi,^ which utilizes Svartholm's^ expansions of Heun functions 
as series of degenerate hypergeometric functions, essentially Jacobi polynomials. First of 
all, it is worth noting that in Erdelyi's method, the root-finding approach is the same as 
that outlined above, because his procedure assumes prior knowledge of the eigenvalues. 
With the eigenvalues already determined, the calculation of the coefl&cients for the Jacobi 
expansion proceeds via a three-term recursion relation similar to that derived by Heun for 
the coefficients of the power-series expansion. Using the familiar result for the quadratic 



hn{x) 



Hn{x) 



(32) 



7-1/2 



with normalization 




(33) 
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normalization integrals of the Jacobi polynomials, it is a simple matter to determine In 
from the coefficients of the Jacobi expansion. While Erdelyi's method is interesting from 
the point of view of functional analysis, it is obvious that his procedure entails much more 
computation than the evaluation of using (29), which we again emphasize requires 
nothing more than information available as a by-product of the search for the eigenvalues. 
In conclusion, we point out that the application of L'Hopital's rule used here to evaluate 
the normalization integrals of Heun functions can also be used to obtain similar results for 
more general Sturm-Liouville problems. 

REFERENCES 

■"^K. Heun, Zer Theorie der Riemann'schen Functionen zweiter Ordnung mit vier 
Verzweigungspunkten, Math. Ann., 33, pp. 161-179 (1889). 

^A. Ronveaiix, ed., Heun's Differential Equations (Oxford University Press, New York, 
NY, 1995). 

^A. Erdelyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi, Higher Transcendental 
Functions, vol. Ill (McGraw-HiU, New York, NY, 1955). 

^A. Erdelyi, Certain expansions of solutions of the Heun equation, Quart. J. Math. 
Oxford Ser., 15, pp. 62-69 (1944). 

^C. G. Lambe and D. R. Ward, Some differential equations and associated integral 
equations. Quart. J. Math. Oxford Ser., 5, pp. 81-97 (1934). 

^N. Svartholm, Die Losung der Fuchsschen Differentialgleichung zweiter Ordnung 
durch hypergeometrische Polynome, Math. Ann., 116, pp. 413-421 (1939). 



12 



